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ABSTRACT 

Angular momentum transport and the formation of rotationally supported structures are major issues in our understanding of protostel- 
lar core formation. Whereas purely hydrodynamical simulations lead to large Keplerian disks, ideal magnetohydrodynamics (MF[D) 
models yield the opposite result, with essentially no disk formation. This stems from the flux-freezing condition in ideal MHD, which 
leads to strong magnetic braking. 

In this paper, we provide a more accurate description of the evolution of the magnetic flux redistribution by including resistive terms 
in the MHD equations. We focus more particularly on the effect of ambipolar diffusion on the properties of the first Larson core 
and its surrounding structure, exploring various initial magnetisations and magnetic held versus rotation axis orientations of a 1 M© 
collapsing prestellar dense core. 

We used the non-ideal magnetohydrodynamics version of the adaptive mesh refinement code RAMSES to carry out these calculations. 
The resistivities required to calculate the ambipolar diffusion terms were computed using a reduced chemical network of charged, 
neutral, and grain species. 

Including ambipolar diffusion leads to the formation of a magnetic diffusion barrier (also known as the decoupling stage) in the 
vicinity of the core, which prevents accumulation of magnetic flux in and around the core and amplification of the held above 0.1 G. 
The mass and radius of the first Larson core, however, remain similar between ideal and non-ideal MHD models. This diffusion 
plateau, preventing further amplification of the held and reorganising the held topology, has crucial consequences for magnetic 
braking processes, allowing the formation of disk structures. Magnetically supported outflows launched in ideal MHD models are 
weakened or even disappear when using non-ideal MHD. In contrast to ideal MHD calculations, misalignment between the initial 
rotation axis and the magnetic held direction does not significantly affect the results for a given magnetisation, showing that the 
physical dissipation processes truly dominate numerical diffusion. 

We demonstrate severe limits of the ideal MHD formalism; it yields unphysical behaviours in the long-term evolution of the system. 
This includes counter-rotation inside the outflow or magnetic tower, interchange instabilities, and flux redistribution triggered by 
numerical diffusion. These effects are not observed in non-ideal MHD. Disks with Keplerian velocity profiles are found to form around 
the protostar in all our non-ideal MHD simulations, with a final mass and size that strongly depend on the initial magnetisation. This 
ranges from a few 10"^ M© and ~ 20-30 au for the most magnetised case (// = 2) to ~ 2 x 10"^ M© and ~ 40-80 au for a lower 
magnetisation (jn = 5). In all cases, these disks remain significantly smaller than disks found in pure hydrodynamical simulations. 
Ambipolar diffusion thus bears a crucial impact on the regulation of magnetic flux and angular momentum transport during the 
collapse of a prestellar core and the formation of the resulting protostellar core-disk system, enabling the formation and growth of 
rotationally supported structures. 
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1. Introduction 

In the past few decades, star formation studies have been strug¬ 
gling to properly describe the mechanism of angular momen¬ 
tum transport, regulated by angular momentum conservation 
and magnetic braking. The formation of rotationally supported 
structures, that is, the protoplanetary disks that are expected to 
give birth to planets and/or binaries, is directly related to the 
amount of angular momentum in protostellar systems. While 
there is plenty of observational evidence for large Keplerian 
disks^ around Class-II and Class-I objects, their presence around 
younger Class-0 objects is still subject to debate (as discussed in 
Maury et al. 2010 or Tobin et al. 2012, 2013, see the review by 
Li et al. 2014a). 

^ See Appendix C for details on our definition of a Keplerian disk. 


The simplified framework of the ideal magnetohydrodynam¬ 
ics (MHD) used in the first studies led to the disappearance of 
the large Keplerian disks that are easily formed in hydrodynami¬ 
cal simulations as a result of the very effective magnetic braking 
created by a strong pile-up of magnetic flux towards the centre 
of the collapsing system (Galli et al. 2006, Allen et al. 2003, 
Price & Bate 2007; Hennebelle & Teyssier 2008; Matsumoto 
& Tomisaka 2004; Hennebelle & Fromang 2008; Commer 9 on 
et al. 2010). In these simulations, disks were found to form 
only for unrealistically weak magnetic field intensities (corre¬ 
sponding to a mass-to-flux ratio more than 10 times the critical 
value derived by Mouschovias & Spitzer 1976). Other conse¬ 
quences of the magnetic flux freezing assumption inherent to the 
ideal MHD approximation include the strong resulting magneti¬ 
sation of protostars compared to the low observed values in stars 
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(Crutcher 2012), the generation of violent interchange instabil¬ 
ities at the protostar-disk interface (Li et al. 2014b), the growth 
of the pseudo-disk (Hennebelle & Fromang 2008), and the dis¬ 
tortion of the upper layer of the pseudo-disk due to reconnection 
and the split monopole (Galli & Shu 1993; Li et al. 2014b). 

In the framework of ideal MHD, there is no possibility to 
regulate the magnetic flux pile-up and its consequences, except 
for the intrinsic numerical resistivity due to both the numerical 
method used to solve the induction equation and the numerical 
grid resolution. To circumvent this problem, the magnetic field 
redistribution needs to be correctly addressed in the framework 
of complete non-ideal MHD. Following the pioneering work of 
Mestel & Spitzer (1956), the study of non-ideal MHD effects 
has been the focus of intense research in the recent years, as il¬ 
lustrated by the studies of Duffin & Pudritz (2008) and Mellon & 
Li (2009) for ambipolar diffusion, or Machida et al. (2006) and 
Machida & Matsumoto (2011) for a generic resistivity. While 
it is still unclear if non-ideal MHD can, by itself, solve the 
problems regarding the formation of disks (Krasnopolsky et al. 
2010), a physical dissipative scale for the magnetic flux defi¬ 
nitely improves the regulation of magnetic flux pile-up in star 
formation studies. 

In this work, we study the effects of ambipolar diffusion in 
the context of star formation, more specifically, of low-mass star 
formation, and highlight the differences compared to ideal MHD 
simulations. The influence of turbulent initial conditions will be 
studied in a forthcoming paper (Paper II). The article is organ¬ 
ised as follows. In Sect. 2 we discuss the framework and nu¬ 
merical setup of the study. In Sect. 3 we focus on the general 
description and properties of the collapsing core until formation 
and early evolution of the first Larson core, for various cases. 
In Sect. 5 we discuss the long-term evolution of the structures 
and highlight the limits of the ideal MHD framework that are 
due to numerical issues. In Sect. 4, we examine the formation 
of rotationally supported structures in non-ideal MHD. Last, we 
summarise our findings in Sect. 7. 

2. General context 


The symbols Vn, p, and pi denote the velocity of neutral parti¬ 
cles, the total density, and the ion density, respectively. In this 
framework, the drag force per unit volume exerted by the neu¬ 
trals on the ions compensates for the Lorentz force and reads 
Fdrag = 7ADP/P(vi - Vn). When considering multiple species (see 
Kunz & Mouschovias 2009; Nakano et al. 2002), Ohm’s law can 
be written as 

Vn XB 
-77a(V X B) 

-,xDj|jix{(VxB)xA}]. (5) 


dB 

dt 


--Vx 


with II II standing for the L2 norm, and the Ohmic, Hall, and 
ambipolar diffusivities are defined as 
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The parallel, perpendicular, and Hall conductivity components 
(cr||, o-±, and (Th, respectively) compose the conductivity tensor 


f 


cr = 


o-H 

0 


-o-H 0 'I 

cr^ 0 

0 cr|| , 


(9) 


from Faraday’s law j = (t(E -r v x B). These resistivities re¬ 
duce to the values displayed in Eq. (2) in a three-fluid descrip¬ 
tion. In this case. 


2.1. Physical framework 

In star formation, the ionisation fraction is low, and quasi¬ 
equilibrium holds between the Lorentz force and the plasma- 
neutrals friction force as a result of collisions. Therefore, we 
can drop the pressure and gravitational forces for charged par¬ 
ticles when writing the equations of motion for each fiuid parti¬ 
cle. In the case of positively charged species of number density 
Hi, velocity vi, atomic number Zi , and collision rate Vij with the 
species j, and using the subscript e for the negatively charged 
species, the equations of motion read: 


1 ZietliiE -1- Vi X B) - Pi i;^=(e,n) n/Vi - \j) = 

\ -ene(E-l-VeXB)-pei;j=(i,„)Ve/Ve-Vy) = 

o o 

(1) 

while the generalized Ohm’s law is written (see e.g. Balbus & 
Terquem 2001) 
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cr\\ = -^- the electrical conductivity. 
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The present work is essentially devoted to the first collapse 
and the formation of the first core. At this stage and on this scale, 
only ambipolar diffusion plays a role. In the following, we there¬ 
fore only focus on this term. Rewriting the above equations with 
? 7 a = ??H = 0 and 


V = v„ -I- 


?7ad 

IIBIP 


[(VxB)xB], 


(13) 


Eq. (5) reduces to 


dfB = V X l^v X Bj. 


(14) 


In contrast to a Laplace operator, as in Ohmic dissipation, 
Eq. (14) does not allow magnetic reconnection, since it corre¬ 
sponds to another flux-freezing condition at a different speed v. 
We stress that this is true only under the one-fluid approximation, 
as pointed out by Tsap et al. (2012). 
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2.2. Magnetic braking 

The main effect of magnetic braking in the context of prestellar 
core collapse is to slow down the rotationally supported struc¬ 
tures that develop during the gravitational collapse of matter with 
non-zero net angular momentum; this occurs because of the con¬ 
servation of angular momentum. One of the consequences is that 
the formation of large and rapidly rotating disks is hampered, as 
found in purely hydrodynamical simulations. The braking arises 
essentially from the magnetic tension that prevents the distor¬ 
tion of the field lines, which are coupled to the flow through the 
ionised species. The angular momentum is transported along the 
held lines at a speed close to the Alfven speed, depending on 
the topology of the magnetic held. The braking efficiency can 
be characterised by deriving an Alembert equation of propaga¬ 
tion for the angular momentum (e.g. Gillis et al. 1974, 1979; 
Mouschovias 1978,1979; Mouschovias & Paleologou 1980), but 
can be estimated by order-of-magnitude calculations (Joos et al. 
2012). In cylindrical coordinates, with the z-axis aligned with the 
mean angular momentum direction vector, the angular momen¬ 
tum flux due to the magnetic held reads (Joos et al. 2012) 




(15) 


It is proportional to the radial and poloidal components of the 
held Br and while scaling with the square of the toroidal 
component Bq, making this latter the main focus of magnetic 
braking studies. 


2.3. Numerical setup 
2.3.1. Methods 

To carry out our study of collapsing magnetised molecular 
cloud cores, we used the adaptive mesh refinement (AMR) code 
RAMSES (Teyssier 2002; Fromang et al. 2006) in its non-ideal 
MHD extension (Masson et al. 2012). RAMSES solves the com¬ 
plete set of MHD equations (self-gravity, Euler’s fluid flow equa¬ 
tions, and the induction equation with non-ideal terms) using the 
constrained transport method, which preserves the divergence- 
free condition for the magnetic held to machine precision. The 
adaptive mesh is extremely well suited to protostellar collapse 
calculations, where many levels of refinement are needed to effi¬ 
ciently describe spatial scales spanning 10"^ - 10^ orders of mag¬ 
nitude in a single simulation. AMR is also a powerful tool for 
fragmentation studies and disk formation in a turbulent medium 
where nested grids are difficult to use. 

The grid refinement criterion is based on the Jeans mass, ensur¬ 
ing the Jeans length is always sampled by at least eight cells. The 
coarse grid has a resolution of 32^, and 11 levels of AMR were 
used, resulting in a maximum resolution of 0.15 au at the finest 
level. Our general set of equations includes the conservation of 
mass, the mean neutral gas dynamics (Euler equation), the induc¬ 
tion equation with the ambipolar diffusion electromotive force, 
self-gravity through the Poisson equation, and the divergence- 
free constraint: 
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^-Vx(vxB)-Vx Ead = 0 , 
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V • B = 0. 

(19) 


Here p is the mean fluid density, v the mean fluid velocity, 
P the thermal pressure of the gas, and Ead = -^ADyfy x 

|(V X B) X jjljjl is the electromagnetic force due to the ambipolar 
diffusion, as derived in Eq. (5). The energy equation is approxi¬ 
mated by a barotropic equation of state (see below). 

Magnetic resistivities are calculated using a reduced chem¬ 
ical network including neutral and charged species, as well as 
dust grains. We followed Kunz & Mouschovias (2009) to com¬ 
pute the relevant charged species abundances including grains 
sizes in a classical MRN distribution (Mathis et al. 1977), which 
were sampled using 50 bins. For an exhaustive description of 
the chemical model used and its application in the context of 
star formation, we refer to Marchand et al. (2015). We computed 
a three-dimensional table of density, temperature, and magnetic 
held dependent resistivities covering the ranges 10‘ -24 < p < 
10“^^ g cm“2, 5 < T < 2000 K, and 10“^ < B < 10^ G, respec¬ 
tively. During the simulations, the resistivities in each grid cell 
are interpolated on-the-fly according to the local state variables, 
which greatly reduces computational cost but implies thermody¬ 
namical equilibrium. 


2.3.2. Initial conditions 


We adopted initial conditions similar to those in Commer 9 on 
et al. (2010), who followed Boss & Bodenheimer (1979). A mag¬ 
netised uniform-density sphere of molecular gas, rotating about 
the z-axis with solid body rotation, is placed in a surrounding 
medium a hundred times less dense with equal pressure. The 
prestellar core mass has a mass of 1 M©, a radius Rq = 2500 au^ 
and a ratio of rotational over gravitational energy of ySrot = b.bl. 
The magnetic held is initially parallel to, and invariant along, 
the rotation z-axis. The held strength is stronger in a cylinder of 
radius Rq (with the dense core at its centre) than in the surround¬ 
ing medium, with B^ir > Ro) = B^iRo) * lOO^/^, where the factor 
of 100 comes from the difference in density between the core 
and the surroundings^. We define a mass-to-flux ratio parameter 
similar to the one defined by Mouschovias & Spitzer (1976) to 
measure the importance of the magnetisation in the core: 





’ 

V 0 /crit 


( 20 ) 


with the critical value ^ (§) (Mouschovias & 

Spitzer 1976). We note that p(r = Ro) is strictly equal to the 
theoretical value for a homogeneous cloud permeated by verti¬ 
cal held lines. 

Even though a Bonnor-Ebert (BE) density profile may better 
fit observations of dense cores (see Andre et al. 2000; Belloche 
et al. 2002) and has an analytical foundation (Li & Shu 1996; 
Hunter 1977), we assumed a magnetised medium permeated by 
straight parallel held lines, and it is hardly possible to end up 
with a BE density profile without bending the held lines during 
the formation of the density enhancement (the dense core). Some 
authors (van Loo et al. 2008, e.g.,) found non-linear density en¬ 
hancements in simulated turbulent molecular sheets via slow¬ 
mode magnetic waves, which left the magnetic held unchanged 

^ The initial condition corresponds to a ratio of the thermal over grav¬ 
itational energies a - 0.25 and a density of 9.4 x 10"^^ g cm"^. 

^ This is chosen to try to reproduce the dragging-in of field lines that 
would occur in the formation of the dense core (see Gillis et al. 1974, for 
example), while also retaining in the simplest manner the divergence- 
free condition for the MHD. 
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in the cores, but it remains unclear whether this process does 
occur in molecular clouds. Our own tests tend to show that the 
initial density profile is not a critical matter because the infalling 
material adopts a BE-like density profile very shortly after the 
beginning of the dense core collapse. 

As the aim of the present paper is to focus on the effect of 
non-ideal MHD effects, notably ambipolar diffusion, on prestel- 
lar core collapse and disk formation, we used a barotropic equa¬ 
tion of state (EOS) to mimic the effect of radiative transfer in¬ 
stead of solving the full set of radiation magnetohydrodynamics 
equations. This was done for the sake of simplicity and to save 
computational cost. The gas pressure is thus related to the den¬ 
sity as 




10“^^ g cm“^ 


4 

3 


( 21 ) 


where Cs is the gas sound speed and nu the number density of 
atoms. After an initial phase of isothermal collapse up to a den¬ 
sity Hu ^ 10“^^ g cm“^, the first hydrostatic Larson core (Larson 
1969) forms when the gas becomes optically thick enough to 
stop radiative cooling, which was until then counter-balancing 
the compressional heating. This adiabatic phase is reproduced 
with the barotropic EOS by using a polytropic index n for a gas 
(which is equal to the ratio of specific heats y) n = |, higher 
than the critical value for stability against gravitational collapse 
^critic = f- The temperature and density of the gas begin to 
rise inside the core as it accretes material from the surrounding 
envelope. The energy sink provided by the dissociation of H 2 
molecules, which occurs around T ~2000 K, triggers a second 
phase of collapse (Larson 1969). We here focus on the proper¬ 
ties of the first Larson core. Our general approach allows us to 
accurately describe the first Larson core and its surroundings and 
the magnetic flux transport without needing to introduce a sink 
particle that would limit the resolution. It is beyond the scope of 
the present paper to include long-term evolution, which leads to 
the formation of the second core. This would imply the accurate 
treatment of jets, high-energy radiation, and remaining non-ideal 
MHD terms and will be addressed in a forthcoming study. 

We have performed eight simulations for which we varied 
the mass-to-fiux ratio p = 2 and 5, the angle between the rota¬ 
tion axis and the magnetic field initial direction (0 and 40 de¬ 
gree) and used ideal (iMHD) or non-ideal (niMHD) magnetohy¬ 
drodynamics (accounting only for ambipolar diffusion). We also 
performed additional simulations with increased and decreased 
resolutions to study the convergence of our results for our fidu¬ 
cial case (p = 5, aligned case). The various run parameters are 
given in Table 1 . 


3. Early evolution 

3.1. Aligned case 

In this section, we study the general properties of the collapsing 
system by comparing the iMHD case and a case with ambipolar 
diffusion in the fiducial aligned case. We first focus on the p = 5 
case, which we describe in detail. We then report the differences 
or similarities with the yu = 2 case. 

3.1.1. First Larson core 

We determined the onset of the first core formation as the 
point when the density in the computational domain reaches 
10“^^ g cm“^. We subsequently defined the first core itself by 


the cells fulfilling the criterion p > 10“^^ g cm“^. The core 
radius, mass, mass-to-fiux ratio, and peak magnetic field value 
were measured 200 years after its formation. The properties of 
the first core are summarised in Table 2. 

The first core is oblate with a radius ranging from 8 to 9 au, 
independently of the initial magnetisation {p = 2 or p = 5). 
While pure hydrodynamical simulations (e.g. Bate 2011; Tomida 
et al. 2010a) and some radiation iMHD (RMHD) calculations 
(Tomida et al. 2010b) found a similar result with a flattened first 
core, other RMHD studies (Commer 9 on et al. 2010) did not find 
flattened cores in the case p = 5. As a result of the additional 
magnetic support, the free-fall (thus core formation) timescale 
is twice longer for p =2 than for /i = 5 . The additional sup¬ 
port also hinders accretion, and the first core is significantly less 
massive for p =2 than for yu = 5 . Similarly, for the stronger 
magnetisation {p = 2), the core mass is lower for niMHD than 
for iMHD. Indeed, the strongest magnetic field inside the core 
remains of the order of 10“^ G in niMHD, while it is ten times 
stronger in iMHD. This effect of ambipolar diffusion on the mag¬ 
netisation during the core formation is illustrated by the differ¬ 
ences in the mass-to-fiux ratio yu at 10 and 100 au between the 
ideal and non-ideal MHD cases. While p is higher than 10 at 
10 au and at < 4 at 100 au in niMHD, it is at most < 5 in iMHD. 


3.1.2. Magnetic field repartition 

Ambipolar diffusion enables neutral particles to overcome mag¬ 
netic field lines and redistribute the magnetic field. We first fo¬ 
cus on the magnetic field intensity and topology to understand 
its consequences for the dynamics of the collapsing core. Lig. 1 
shows the magnetic field repartition as a function of gas density 
for the iMHD and niMHD runs. The flux-freezing behaviour, 
B oc with ^ ~ 3/2, is obvious for iMHD (red) with a 
highest magnetic field value of ~1 G. The niMHD simulation 
(blue) shows an identical initial evolution for densities below 
~ 10“^"^ g cm“^. After this point, ambipolar diffusion starts to 
dominate the later evolution of the core, and a very well defined 
diffusion plateau forms that is usually referred to as the decou¬ 
pling stage in star formation (Desch & Mouschovias 2001). In 
both the iMHD and niMHD calculations, the generation of out¬ 
flows can be identified by the high values of the magnetic field 
at high density, compared with the expected perfect flux-freezing 
result. 

As seen in Lig. 1, a scaling relation B oc seems to repre¬ 
sent the evolution of the magnetic field during the collapse better 
than the traditional B oc ^Jp relation (see the differences in the 
scaling for the different components of the field in Hennebelle & 
Lromang 2008). Lurther details about this scaling relation and 
the analytical derivation of an estimate of the field saturation 
value are given in Appendix A. 

Lig. 2 illustrates the comparison between p = 2 and p = 
5 niMHD. Both simulations show similar initial evolutions 
(Lig. 2a) when flux-freezing is still dominant (p = 2 naturally 
shows a higher magnetic field intensity at a given density), but 
have slightly different saturation values. At later stages (Lig. 2b) 
a plateau forms in both low magnetisation and high magnetisa¬ 
tion runs with ^plateau - 0.1 G. 

This plateau arises when ambipolar diffusion becomes ef¬ 
ficient enough to overcome the dynamical effects. The non¬ 
linearity of the effective diffusion coefficient, 77 ad ex¬ 

plains the self-regulation and the formation of the diffusion 
plateau. Indeed, for any non-potential (V x B 9 ^ 0) field con¬ 
figuration, an increase in the magnetic field leads to a rise of 
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Table 1. Summary of the initial conditions for the different simulations 


Alignment and MHD type 

niMHD Aligned 

iMHD Aligned 

niMHD Misaligned 

iMHD Misaligned 

magnetisation (p) 

2 

5 

2 

5 

2 

5 

2 

5 

Ambipolar diffusion 

yes 

yes 

no 

no 

yes 

yes 

no 

no 

Angle 

0 

0 

0 

0 

40 

40 

40 

40 

Rotational support (jS^ot) 

0.02 

0.02 

0.02 

0.02 

0.02 

0.02 

0.02 

0.02 

Mass of the core (in M©) 

1 

1 

1 

1 

1 

1 

1 

1 


Table 2. Summary of the main properties of the first core for fiducial cases with solid-body rotation 


Alignment and MHD type 
magnetisation 

niMHD Aligned 
p = 2 p = 5 

iMHD Aligned 
p = 2 p = 5 

niMHD Misaligned 
p = 2 p = 5 

iMHD Misaligned 
p = 2 p = 5 

Formation time (xlO^ years) 

49.2 

24.3 

50.5 

24.3 

52.2 

24.6 

54.6 

24.6 

Radius (au) 

8-9 

8-9 

9-10 

7-9 

8-9 

8-9 

9-10 

8-9 

Mass (Mq) 

0.016 

0.036 

0.025 

0.03 

0.016 

0.058 

0.024 

0.058 

p(r = 10 au) 

14.1 

21.2 

3.5 

3.9 

14.1 

24.8 

3.4 

3.9 

p(r =100 au) 

2.8 

4.2 

1.8 

4.6 

2.8 

3.5 

1.6 

2.8 

Magnetic field strength B (G) 

0.09 

0.1 

0.8 

0.98 

0.09 

0.21 

0.77 

4.4 



-18 -16 -14 -12 -10 -8 

log(p) (g cm‘^) 


Fig. 1. Magnetic field distribution as a function of gas density for 
fjL - 5 aligned, using the iMHD (red) and the niMHD (blue) formalisms, 
200 years after first core formation. The solid black line shows the scal¬ 
ing of the magnetic field B oc and the dashed black line the scaling 
B oc . 


the local diffusion coefficient and thus in turn to a decrease in 
field intensity. The plateau first forms when the dynamical evo¬ 
lution of the magnetic field, v x B, no longer dominates the dif¬ 
fusion term, t/adV x B x jjljj x jj|jj (with || || standing for the 
the L2 norm in this paper). Assuming flux-freezing holds during 
the first isothermal phase of the collapse (i.e. B oc we can 
estimate the saturation value for the magnetic field 


^sat — I ^TAD \l 2jiQ 


2-^ 






Pq 


( 22 ) 


where 7 ad is the drag coefficient (Eq. 3), C the ionisation frac¬ 
tion such that Pi = C G the gravitational constant, 
and po the typical initial conditions for the magnetic field and 
gas density, and ^ the power index for the proportionality law 
B(p) oc p^/^. The initial thermal support and mass-to-flux ra¬ 
tio are linked to the values of Bq and po. Further details on the 
derivation of ^sat can be found in Appendix A. The values we ob¬ 
tain for ^sat with ^ = I for various values of thermal support a, 
that is, the ratio of thermal over gravitational energy, are shown 
Fig. 3 (thick black solid line). For ^ which corresponds to a 


homogeneous sphere permeated by a uniform magnetic field, we 
And that the saturation value does not depend on a. In this case, 
the saturation value estimate for p =5 is ^sat = 0-13 G, very 
close to the numerical value 5 < 0.1 G (see Fig. 1). According to 
this simple order-of-magnitude calculation, the saturation value 
should even be lower in more magnetised cores, with p < 5. For 
example, the p = 2 case yields ^sat = 8 x 10“^ G. In our numer¬ 
ical simulation, we do not observe such a strong diminution, but 
the saturation value is lower for p = 2 than for p = 5 (see Fig. 2). 
The discrepancy between the analytical prediction and the nu¬ 
merical value for ^sat in the strongly magnetised case originates 
from the fact that in the latter case the magnetic field distribu¬ 
tion departs appreciably from a simple power-law parametrisa- 
tion B oc p^/^ (see Fig. 2). Using a slightly different power index 
that better fits the efiTective B(p) repartition in the p =2 case, 
for instance ^ = 1.73, as shown on the inset in Fig. 3, we obtain 
^sat = 0-06 (coloured lines in Fig. 3), which is closer to the value 
obtained in the simulation. 

Fig. 4 shows snapshots of the mass-to-flux ratio p, as de¬ 
fined Eq. (20), as a function of radius at different times in the 
simulations. At large radii (r > 1000 au), the magnetisation 
is similar for iMHD (red) and niMHD (black), confirming the 
fact that at these scales the AD timescale is orders of magni¬ 
tude longer than the dynamical timescale. There is a kink in the 
AD case at a few tens of au (highlighted by the green square in 
the figure) that slowly propagates outwards (compare the blue 
dashed and solid curves); this corresponds to the efficient diffu¬ 
sion of the field at these densities. 


3.1.3. Outflows 

The piling-up of the toroidal component of the magnetic field 
during the collapse of the dense core ultimately leads to the for¬ 
mation of a growing vertical structure, called magnetic tower. 
We note that we use the generic term magnetic tower to describe 
any outflowing structure. In reality, there are several ways to 
launch outflows in a magnetised environment, as studied initially 
by Lynden-Bell & Boily (1994) and in contemporary studies by 
Hennebelle & Fromang (2008); Ciardi & Hennebelle (2010). 
This outflow is launched shortly (< 1 kyear) after the forma¬ 
tion of the first core. Fig. 5 (top row) displays the azimuthally 
averaged density and velocity fields (panels a and c) and Alfven 
speed with the magnetic field direction (panels b and d) for the 
iMHD and niMHD p = 5 runs. The simulations are compared 
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log(p) (g cm-'^ 

(a) Early time, when flux freezing is still relevant. 



log(p) (g cm-'^ 

(b) After the decoupling stage. 


Fig. 2. B(p) for // = 5 (red) and p = 2 (blue). 


500 years after the formation of the first Larson core, but the 
same conclusions hold during all the time evolution. In all pan¬ 
els, the z-axis is aligned with the direction of the average angular 
momentum vector computed in a sphere of radius 100 au, cen¬ 
tred around the densest cell in the simulation. The velocity vec¬ 
tors are overlaid on the density maps, while the direction of the 
magnetic field in the (r, z) plane is superimposed on the Alfven 
speed colour maps. The colour coding for the magnetic field il¬ 
lustrates the intensity of the toroidal component 5^. 

The magnetically driven outfiow is reinforced in the iMHD 
case by the more important magnetic field pile-up. This accu¬ 
mulation stems both from the radial bending of the field lines, as 
clearly seen in Fig. 5b, and from the increasing toroidal field (of 
prime importance for magnetic braking, as shown by Eq. 15) that 
is generated by the rotation, as clearly seen in Fig. 5b and d, es¬ 
pecially close to the first Larson core where the toroidal support 
is significantly stronger in (b). While the gas inside the magnetic 
tower is slightly less dense in the AD case, the lower magnetic 
field intensity (compared to the iMHD case) yields an overall 
Alfven velocity that is lower by one order of magnitude. This 
in turn explains the weaker magnetic braking, since the angular 
momentum is carried away by slower Alfven waves^. 


In a simple representation of bent field lines, it is possible to derive 
the angular momentum transport equation. Angular momentum propa- 



Fig. 3. Saturation value as a function of p for various values of thermal 
support a. The points corresponding io p - 2 (lower values) and p - 5 
(upper values) are marked in the figures. 



Fig. 4. p - 5, aligned case. Mass-to-flux ratio (p) after first core for¬ 
mation (at ~ 24 kyears) for the niMHD (black) and IMHD (red) cases 
(solid lines). We also display the initial condition (dotted blue line) and 
two later outputs for the niMHD case (dashed and dotted solid black 
lines). The decoupling stage region due to magnetic field pile-up is em¬ 
phasised with a green square. 


In the iMHD simulation, the velocity field in the tower (in the 
r-z plane) is vertical and the gas mainly fiows out, the launching 
mechanism taking root in the strong toroidal field and differen¬ 
tial rotation close to the Larson core. The field is dominantly 
vertical, and neutral matter follows the field lines. When AD is 
included, the growth of the magnetic tower is still magnetically 
regulated, but the growth is slower and the velocity field in the 
tower is almost null, or slightly directed towards the protostar. A 
detailed analysis shows that there is no real outflowing gas mo¬ 
tion of gas per se, but that the magnetic tower structure or growth 
is supported by magnetic pressure. We also note that close to the 
core (r < 50 au), while the field lines are pinched in the iMHD 
case (split monopole topology), field re-distribution has operated 
in the niMHD run because of the ambipolar diffusion (the mag¬ 
netic field vectors are much more vertical). 


gates along the field lines and follows a wave equation with a velocity 
defined by the local Alfven speed and the topology of the field. For more 
details see Masson (2013) or the original article by Gillis et al. (1979). 
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Fig. 5. Top row: jd - 5, aligned case, snapshots at a time t - 1.02 tfree-faii (500 years after first core formation). Panel (a) shows a map of the density 
in the ideal MHD simulation. The black solid lines are logarithmically spaced density contours. The arrows represent the velocity field where 
blue to green corresponds to the increasing magnitude of the velocity in the (r, z) plane. Panel (b) is a map of the Alfven speed (blue), with the 
magnetic field direction in the (r, z) plane indicated by the segments. Yellow to red in these segments corresponds to increasing relative magnitude 
of the toroidal component of the field compared to the total field. Panels (c) and (d) are the same as (a) and (b), respectively, but for the niMHD 
case. Bottom row: same as the top row for yu = 2, and at the same time t - 1.01 tfree-faii (500 years after first core formation). Every quantity is 
azimuthally averaged. 


li = 2 : In a more magnetised case, the picture is very different. 
The bottom row of Fig. 5 shows the structure of the outflow for 
jd = 2 . Magnetic braking is much more effective (the magnetic 
field is overall stronger, yielding a stronger magnetic braking, as 
seen Sect. 2.2), causing the field topology to come closer to a 
split monopole configuration at large scale, with strong pinching 
of the lines in the equatorial plane. The enhanced braking weak¬ 
ens the outflow-launching mechanism because the pile-up of the 
toroidal field is less effective. For iMHD, there is still an outflow¬ 
ing feature, but its velocity field is either null in the r-z plane or 
even falls back onto the core. The boundary of the outflow re¬ 
gion is less well defined, but is characterised by a discontinuity 
in the velocity field. In the AD case, the braking is stronger for 
jd = 2 than for jd = 5 at the early stages of the collapse (see 
Sect. 3.1.4), which results in a weaker magnetic field (the dif¬ 
fusion plateau is lower) and reduces the pile-up of toroidal field 
close to the core. This weakens the launching process to the point 
that no outflow is produced in this case. The Alfven speed close 
to the core is very similar to the weak field jd = 5 run (Fig. 5d), 


as expected from the field magnitude distribution as a function 
of density (Fig. 2), where B(p > 10“^^ g cm“^) distributions are 
almost identical for the strongly and weakly magnetised cases. 


3.1.4. Regions of active ambipolar action 

Figure 5 also shows that for R >100 au the iMHD and AD 
runs look similar. This resemblance arises from the strong de¬ 
pendence of ambipolar diffusion efficiency upon density and the 
sharp transition between the flux-freezing and AD dominated 
regimes. The initial isothermal phase of the collapse thus re¬ 
mains very similar in both cases. To examine the effect of flux 
dissipation in detail, we now focus on the regions where ambipo¬ 
lar diffusion dominates. 

These regions are highlighted in Fig. 6a, which shows maps 
of Am, an adimensional number that characterises the efficiency 
of the diffusion process compared to the dynamical ones, defined 
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Fig. 6. Four main panels show the // = 5 aligned (6a), p - 5 misaligned (6b), p - 2 aligned (6c), and p - 2 misaligned (6d) cases 500 years 
after first core formation. In each main panel, subpanel (a) shows a colour map of the azimuthal average of Am = ^. Green to blue means that 
dynamical processes (collapse, Keplerian disk) dominate, while red to yellow means that ambipolar diffusion dominates. The grey segments show 
the direction of the magnetic field, and black solid lines are isodensity contours. Subpanels (b) and (c) display the magnetic field magnitude B and 
the value of Am as a function of density. Red cells in the Am scatter plot are cells where B(p) > 0.08 G, while grey cells have B(p) < 0.08 G. 


as 

Am = — . (23) 

^AD 

The typical length-scale L is taken as the distance to the proto¬ 
star and the velocity V as the local velocity along the field lines, 
V = (v • B)/||B||. Ambipolar diffusion essentially plays no role 
in regions where Am > 1, which is the case for r > 100 au, 
where niMHD and iMHD calculations produce identical struc¬ 
tures. The region of the outflow, in contrast, is a region of very 
active ambipolar diffusion (Am < 1) because of a lower den¬ 


sity that corresponds to a higher resistivity t/ad (see Eq. 2). Dis¬ 
sipative effects are also very strong around the mid-plane for 
r < 50 au. This reduces magnetic braking by relaxing the split 
monopole configuration close to the dense core. The right pan¬ 
els in Fig. 6 show the magnetic field distribution as a function 
of density, along with the Am number (right axis). All cells 
with B > 0.08 G (horizontal dot-dashed line) have Am values 
coloured in red (compared to non-flagged cells, which are grey). 
We note that almost every cell belonging to the magnetic field 
plateau is indeed dominated by ambipolar diffusion effects rather 
than dynamical effects (Am < 1). 
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li = 2 : The same diagram for the more magnetised case is 
shown in Fig. 6b. Dynamical motions are less dominant because 
of the additional magnetic support. As a consequence, the active 
ambipolar diffusion region is stretched into an hourglass shape 
with a pinched equatorial waist in the mid-plane for r < 50 au, 
of about the same size as for yu = 5 . The pinching of the field 
lines is also weaker in this central region than for iMHD. How¬ 
ever, compared to yu = 5, the split monopole configuration re¬ 
mains, with a more radial field in the midplane, especially in the 
domain 50 < r < 150 au. This eventually leads to increased 
magnetic braking, which significantly hampers the formation of 
large Keplerian disks. Cells belonging to the diffusion plateau 
are, as previously, dominated by ambipolar diffusion effects and 
characterised by Am < 1. 


3.2. Misaligned case 

The effects of tilting the rotation axis of the molecular cloud 
dense core with respect to the orientation of the magnetic field 
has consequences on the properties of the first Larson core and 
the accretion disks that can form around it (see Joos et al. 2012). 
We have repeated this analysis in MHD calculations including 
ambipolar diffusion. 


3.2.1. First Larson core 

The properties of the first Larson core in the misaligned case 
are summarised in Table 2. As in Sect. 3.1.1, the listed quan¬ 
tities were measured 200 years after peak density reached 
10“^^ g cm“^. 

In the misaligned case, the regulation of the magnetic flux 
occurs as in the aligned case, yielding similar properties for the 
mass-to-flux ratio repartition and the value of the strongest mag¬ 
netic field. For the collapse phase, the additional magnetic ten¬ 
sion leads to a slightly longer free-fall time. For the weakly mag¬ 
netised yu = 5 , the mass of the core is increased in the misaligned 
case. Indeed, the mass-to-flux ratio reaches the same value, but 
the magnetic field is twice stronger than in the aligned case, 
yielding a core mass almost twice higher. 

Therefore, misalignment of the magnetic field and rotation 
axis does not change the properties of the first core significantly. 
It affects the formation and properties of rotationally supported 
structures, however, because of the weaker magnetic breaking, 
as we examine in Sect. 4.2 


3.2.2. Magnetic field repartition 



-18 -16 -14 -12 -10 -8 

log(p) (g cm'^) 

(a) yu = 5 


0 - iMHD 



-5 - ^^^^^ - 

-18 -16 -14 -12 -10 -8 

log(p) (g cm'^) 

(h)fi = 2 


Fig. 7. Same as Fig. 1 for the misaligned case. 


3.2.3. Outflows 

We consider only niMHD in this paragraph. Because of the mis¬ 
alignment angle, outflows, when present, exhibit a wider open¬ 
ing angle in the misaligned case than in the aligned one, and 
their growth is hindered. Misalignment yields less pile-up of the 
toroidal field, which weakens the outflow-launching mechanism. 
As for the aligned case, we do not see any outflowing motions in 
the yu = 2 calculations. 


Magnetic field repartition as a function of density is qualita¬ 
tively unchanged compared to the aligned case, as seen in Fig. 7. 
The same non-linear self regulation process yields a diffusion 
plateau. The region above and below the core is strongly de¬ 
pleted, leading to a highly magnetised and low-density region 
(p < 10“^^ g cm“^, B > 10“^ G) that corresponds to the polar 
cavity mentioned in Tomida et al. (2015). It is more extended 
in the more magnetised case (see Fig. 7a compared to Fig. 7b), 
as in the aligned configuration (see Figs 6c and 6d). In iMHD, 
flux freezing still yields high peak magnetic field values, but the 
mixing due to the misalignment increases the numerical recon¬ 
nection. As a result, a plateau somewhat similar to the niMHD 
case starts to develop (especially in the yu = 2 case). We will 
discuss this in further detail in a forthcoming paper (Paper II), 
which will include the effect of turbulence. 


We conclude that the presence or absence of outflowing gas 
is closely linked to conservation of angular momentum and to 
the regulation of magnetic field accumulation (and thus magnetic 
braking). Low magnetisation (fi = 5) allows for enough toroidal 
field accumulation to grow a magnetically and density-enhanced 
tower, while higher magnetisation (p = 2) produces too effective 
braking, preventing the formation of magnetically launched and 
supported outflows. 


3.2.4. Region of active ambipolar action 

Fig. 6c shows the efficiency of the ambipolar diffusion along 
with density contours and magnetic field orientation for yu = 5. 
The results are similar to the aligned case (cf. Fig. 6a), with a 
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larger pseudo-disk^. At these mid- to large scales (> 50 au), am- 
bipolar diffusion has little influence (Am » 1). The polar cavi¬ 
ties above and below the core are still present and well defined, 
and the region of ambipolar action in the mid-plane leads to the 
formation of a rotationally supported structure. Again, the mag¬ 
netic plateau that forms aX B < 0.1 G corresponds to regions 
where Am < 1, dominated by ambipolar diffusion. 

The = 2 misaligned case is almost indistinguishable from 
the aligned one, as can be seen when comparing Figs 6b and 6d. 
This suggests that in more magnetised cases the final configura¬ 
tion of the protostellar system results essentially from a satura¬ 
tion process (self-regulation from the ambipolar diffusion) and is 
independent of the alignment or misalignment. Another way to 
formulate this result is that all the relevant regions, including the 
diffusion plateau and most of the simulated box, are dominated 
by ambipolar diffusion (Am < 1). 

4. Long-term evolution in non-ideai MHD: disk 
formation 

In this section, we examine the rotationally supported structures 
that can form around the first Larson core. We focus on niMHD 
calculations; iMHD calculations are discussed in Sect 5. 

4.1. Aligned case 

Conservation of angular momentum during the collapse yields 
a very high specific angular momentum close to the protostar, 
which eventually leads to the formation of rotationally domi¬ 
nated structures with a strong toroidal velocity component and 
in some cases to genuine disks with Keplerian velocity profiles. 
To define a "disk" in the simulations, we used the criteria defined 
in Joos et al. (2013). Structures in which the magnetic pressure 
either dominates or is at least not negligible, however, may sat¬ 
isfy these criteria while being very different from the flat rota¬ 
tionally supported structures characteristic of Keplerian disks. It 
is important to clearly identify and distinguish these two types 
of structures, since the existence of flat disks around Class-0 ob¬ 
jects remains a subject of heated debate, and this issue is ad¬ 
dressed below. Following Joos et al. (2013), a piece of fluid be¬ 
longs to the disk if it fulfils all the following constraints (using 
/ = 2 ): 

- it is close to hydrostatic equilibrium: ve > fv^ ; 

- it has strong rotational support: vq > fvy ; 

“ its rotational support is stronger than the thermal support: 

- it has a high density: p > 3.8 x 10“^^ g cm“^. 

In our fiducial aligned case with p = 5, we do observe the 
formation of a disk according to the above criteria that strongly 
resembles a Keplerian disk. Its characteristic properties are given 
in Table 3 and the results of the simulations are shown in Fig. 8 
(top row). The disk is well defined in density, beginning just at 
the edge of the core and exhibiting sharp edges at its periphery. 
The breaking of symmetry observed in the top view (Fig. 8b) 
is due to the evolution over several orbital periods at the core 
radius, with small numerical errors leading to a bar-like instabil¬ 
ity that evolves into two spiral arms. The disk then grows to a 
radius of r ~ 80 au and is flat. The velocity profile in the disk 

^ The pseudo-disk is the overdensity region resulting from the loading 
of pinched field lines in a collapsing core that resembles a disk. See 
Galli & Shu (1993) for the first study of the pseudo-disk. 


Table 3. Properties of the disks that are formed in the non-ideal MHD 
simulations. 


Alignment 

magnetisation 

Aligned 

p = 2 p = 5 

Misaligned 
p-2 p- 5 

Disk age (k years) 

5.1 

2.1 

7.6 

2.5 

Disk (aspect ratio) 

core (2) 

spiral (8) 

spiral (5) 

warped (5) 

Outer radius (au) 

15 

80 

30 

45 

Mflrstcore (Mq) 

0.17 

0.23 

0.17 

0.27 

Mfiisk (Mq) 

0.029 

0.14 

0.027 

0.19 

yS(inner radius) 

> 10^ 

> 10^ 

> 10^ 

> 10^ 

/5(outer radius) 

< 1 

~ 1 

~ 1 

< 1 



Fig. 9. p = 5, aligned case. Toroidal (v^) velocity as a function of radius 
for disk cells. The pink dashed line is the fit (in the log/log space) of 
the blue points by a power law a:^. The grey dashed line is the same fit 
but by enforcing b = -1/2. The vertical dotted line at 7 au is the inner 
radius for the fit. 


is displayed Fig. 9 and is very close to a Keplerian one, with 
an estimated mass for the central object of 0.17 M© compared 
to 0.23 M© (Table 3) (see Appendix C for an important remark 
about assuming a Keplerian profile to estimate the mass of the 
central object). 

The evolution of the disk for our two values of magnetisation 
is displayed in Fig. 10 (solid lines for p = 5, dashed lines for 
p = 2). In niMHD, the pinching of the field lines the midplane is 
greatly reduced and the radial component of the field is almost 
null. In this case, there is almost no magnetic flux accreted onto 
the central object in the radial direction. 

Figures 8c and 8d show the plasma p, defined as the ratio 
of thermal over magnetic pressure, p = F/Fmag. in the disk and 
its surroundings. We note that everywhere inside the disk and 
the core, » 1, meaning that the thermal pressure dominates 
magnetic pressure in these regions. 

These results are of prime importance both from a physical 
and numerical point of view. Physically, they bear major conse¬ 
quences on our understanding of fragmentation and angular mo¬ 
mentum transport in disks. In ideal MHD calculations, magnetic 
fields have been shown to inhibit fragmentation (Hennebelle & 
Teyssier (2008) and Commer 9 on et al. (2010). The fact that 
in more realistic non-ideal MHD calculations yS » 1 suggests 
that magnetic fields have less impact than expected according to 
iMHD results and that fragmentation may be facilitated. A sec¬ 
ond important consequence of the present calculations, of the 
characterisation of plasma p, and of the detailed topology of the 
field is the major role of these quantities in determining viscous 
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Fig. 8. Visualisation of the disk structures. Top row: Panels (a) and (b) show side and top views, respectively, of the disk density (orange), with the 
velocity vectors superimposed for yu = 5 . Panels (c) and (d) are the same side and top views, but showing the plasma (5 parameter (colour map), 
over which we plot the magnetic field direction for the side view (c) and the velocity vectors for the top view (d). Bottom row: Same as for the top 
row, but for the yu = 2 simulation. Each row has a different spatial scale. 



Fig. 10. Disk mass evolution. Solid lines: jd - 5\ dashed lines: ju = 2. 
Black: aligned case. Red: misaligned case, t = 0 corresponds to the first 
core formation to allow direct comparison between cases. 


transport in disks, in particular for the magneto-rotational insta¬ 
bility (Balbus & Terquem 2001). A better knowledge of these 
quantities will enable us to define more accurate initial condi¬ 
tions in such studies (see e.g. Lesur et al. 2014). 

Form the numerical point of view, the present studies are 
very important for the use of sink particles in simulations. Cur¬ 
rent sink particle models do not treat magnetic flux transport at 
the sink-gas interface in a consistent way. The flux is not ac¬ 


creted, which creates an effective barrier of infinite diffusion and 
spurious numerically driven interchange instabilities. As men¬ 
tioned above, in the AD case, the radial component of the mag¬ 
netic field is almost non-existent, yielding a configuration of zero 
net radial flux, allowing the accurate characterisation of the ac¬ 
creted (or in this case not accreted) magnetic flux onto the sink 
particle. 


Id = 2 : In the more magnetised case, a disk-like structure 
around the first core is still observed, according to the above cri¬ 
teria, but the disk mass is an order of magnitude lower than in 
the less magnetised jd = 5 case. Figures 8e and 8f represent the 
disk cells in an edge-on and in a top view. As for the less mag¬ 
netised case, the disk inner radius coincides with the outer layers 
of the core, and the disk slowly grows as matter is accreted. The 
final radius is about 20 au with a mass of 0.03 M©. The aspect 
ratio remains close to unity, and the structure exhibits the same 
characteristic properties as in the previous case: it has a close-to- 
Keplerian velocity field with a high plasma p in the rotationally 
dominated structure, while in contrast, the regions outside the 
disk are magnetically dominated. No spiral arms have developed 
in this case, even at the end of the simulation, but steady state 
has not been reached yet as the disk is still slowly accreting. 


As an important global remark, we note that since values of 
yS » 1 in the disk are found in all our simulations in niMHD for 
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any initial condition, a yS-based criterion could thus be used as a 
robust disk criterion in niMHD simulations^. 

4.2. Misaligned case 

In the misaligned case, a disk can form more easily and is more 
massive, as already noted in previous studies (Joos et al. 2012; Li 
et al. 2013). Indeed, misalignment yields a thicker pseudo-disk, 
which reduces the magnetic braking. Fig. 11 (top row) shows the 
disk structure at the end of the simulation. It exhibits a flat inner 
morphology, with a Keplerian velocity profile and high plasma p 
values. In this case we note large accretion spiral arms (with sig¬ 
nificant infalling motions) above and below the disk plane, with 
P These spiral arms, however, do not represent a significant 
part of the disk mass. We note that using a/5-based disk selection 
criterion would exclude these arms from the disk (see Fig. lid 
compared to 11b). 

For p = 2 (Figs, lle-h), the early evolution is similar to the 
aligned case, the disk is at first hardly distinguishable from the 
core (outer radius is about 10 to 15 au) and then grows by accre¬ 
tion. Significant differences compared to the aligned case, how¬ 
ever, occur at later stages, as seen in the figures. In this case, the 
disk is massive enough to trigger the formation of spiral arms, 
with a mean radius of 20 to 50 au. As for the previous cases, the 
disk is very well defined by the region of high p plasma. 

4.3. Emerging consistent picture 

Fig. 12 shows a direct visual comparison of disks in each simula¬ 
tion. The disk formation and evolution during the collapse in the 
aligned and misaligned cases is presented in Fig. 10. In the mis¬ 
aligned cases, the disk forms at the same time as or even before 
the first core forms. This confirms, as discussed above, that the 
disk emerges from the outer layers of a distorted first core em¬ 
bryo. After about a hundred years, the disk is well defined, with 
a mass > 10“^ Mq in all cases. At this stage, the disk growth 
histories for the aligned and misaligned cases are essentially in¬ 
distinguishable. 

Fig. 13 displays the angular momentum evolution in spatially 
restricted regions, separating the disk and outflow components in 
the system. The disk plane includes every cell fulfilling the disk 
criteria inside a sphere of radius 100 au. Conversely, the outflow 
component corresponds to every other cell in a sphere of radius 
100 au. The p = 5 and p = 2 cases exhibit the same pattern. 
The mass-averaged angular momentum is slightly larger in the 
outflow region in the aligned case (dashed blue line) than in the 
misaligned one (dashed red line) during the first thousand years 
after the core formation. In the disk plane, the definition of the 
disk and the presence of spiral arms introduce variability, but the 
trend remains similar for the aligned and misaligned configura¬ 
tions. The misaligned case in iMHD releases the strong gradients 
in the pseudo-disk (strongly pinched magnetic held lines) found 
in the aligned case and thus enables larger disks to form. Non¬ 
ideal MHD produces the same smoothing of gradients, yielding 
eventually similar structures in the aligned and misaligned cases. 

This similarity between aligned and misaligned simulations 
in niMHD shows that the formation and evolution of a disk is 
independent of the initial misalignment. Ambipolar diffusion, 
when it is efficient, regulates the angular momentum transport 
and the magnetisation during the collapse. A stronger magnetic 

^ This criterion encompasses the core itself, which needs to be re¬ 
moved from the disk by considering the thermal to kinetic energy ratio, 
as done in Joos et al. (2013). 



0 50 100 150 

Radius (au) 


Fig. 12. Representation of disk sizes and shapes for every niMHD sim¬ 
ulation. The density colour map and contours are the same as in Figs. 8 
and 11. 



Fig. 13. Angular momentum evolution. Blue and red: p = 5; grey and 
black: p - 2. Solid lines: angular momentum in the disk plane; dashed 
lines: angular momentum in the outflow (see text). ^ = 0 corresponds to 
the formation of the first core. 
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Fig. 11. Same as Fig. 8, but for the misaligned configuration. Each row has a different spatial scale. 


braking will operate in the aligned case, leading to more angu¬ 
lar momentum being evacuated in the vertical direction and ulti¬ 
mately yielding the same disk mass and final angular momentum 
once a physical equilibrium has been reached. The larger the am¬ 
bipolar diffusion, the more efficient the regulation mechanism, 
leading eventually to similar structures. This is well illustrated 
by the even more pronounced similarity between the aligned and 
misaligned cases for jj = 2. 

5. Limits of ideai MHD 

In this section, we highlight several limits of the ideal MHD ap¬ 
proximation in the context of prestellar magnetised collapse and 
disk formation, justifying in passing the fact that we did not dis¬ 
cuss iMHD simulations in Sect. 4. 

5.1. Counter-rotation 

In ideal MHD, the increase in the toroidal field component is 
due to the rotation of the gas around the forming protostar and 
is only hindered by numerical reconnection. The resulting mag¬ 
netic braking can be efficient enough to completely stall the rota¬ 
tion of the core. During the time the information propagates out¬ 
wards, the core can start to counter-rotate, creating a new brak¬ 
ing, until eventually co-rotation with its surrounding is reached 
again. This is illustrated in Fig. 14, where we slightly increased 
the initial rotation in our fiducial case to better illustrate our pur¬ 
pose: fi = 5, a = 0.25,ySrot. = 0 - 03 . Time is evolving from top 
to bottom. To trace the rotation, we calculated the angular mo¬ 
mentum for each cell with respect to the mean direction of rota¬ 
tion in a sphere of radius 200 au. The first and second columns 
show the negative and positive toroidal velocity {vq) maps, re¬ 


spectively, in a iMHD simulation. The third column displays the 
positive toroidal velocity for niMHD Strong counter-rotation 
occurs in the iMHD simulation, starting near the core (second 
row), and propagates into the outflow (third row). There is no 
sign of counter-rotation in the niMHD run. The mechanism con¬ 
tinues and generates zones that extend to very large radii with 
alternate positive and negative toroidal velocities, producing a 
butterfly-shaped outflow (bottom row). This has important con¬ 
sequences on the expansion of the outflow, which is narrower for 
iMHD (this is most visible in the third row) because the counter¬ 
rotation hampers the pile-up of the toroidal component of the 
magnetic field. 

We carried out similar studies in the misaligned case and 
found out that counter-rotation still develops after formation of 
the first core, although with a more limited spatial extent than 
in the aligned case. We also examined lower initial magnetisa¬ 
tions and found out that counter-rotation is greatly reduced when 
yU > 5. 

5.2. Flux redistribution 

Fig. 15 illustrates the evolution of the disk mass and central mag¬ 
netic field for iMHD . At ^ > 26.5 k years, the field intensity 
decreases by almost an order of magnitude, producing a drastic 
increase of the disk mass. This stems from the flux accumulation 
at the centre of the collapsing system because the magnetic flux 
is frozen with the flow, which produces a major interchange in¬ 
stability. After this event, the core environment is weakly magne¬ 
tised and quite disorganised (looking similar to turbulent runs), 

^ We do not display the counter-rotating < 0 for AD since there are 
almost no counter-rotating cells except very close to the rotation axis at 
the very end of the simulation. 
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IMHD IMHD AD 



Fig. 14. p - 5, aligned case. Four snapshots (time increases from top 
to bottom) from an ideal MHD simulation (two left columns) and the 
ambipolar diffusion case (right column) with the magnetic field initially 
aligned with the rotation axis. For the iMFID simulation, the angular 
momentum is plotted in the left or right panel, depending on the sign of 
the azimuthal speed vg in cylindrical coordinates. 


enabling the formation of a massive rotationally supported disk. 
The evolutionary sequence is portrayed in Fig. 16 (top), where 
the times corresponding to snapshots 1 to 4 in Fig. 15 are indi¬ 
cated in Fig. 15. The thin disk-like structure in 1 and 2, with a 
mass Mdisk ^ 5 x 10“^ Mq, evolves into a thick structure (3 and 
4) with strong magnetic support, although it does not fulfil the 
disk criteria we defined above. The violent release of magnetic 
flux produced by the interchange instability between epochs 2 
and 3 breaks the top-down symmetry of the system and displaces 
the core by several tens of au (this is not visible in Fig. 15 since 
we have integrated the maps by always using the densest cells as 
the origin). Detailed maps of the different magnetic components 
and the plasma p are displayed in Fig. 16 (second and third rows) 
for snapshots 1 and 4. The symmetry breaking, as well as the 
toroidal support, are apparent in the figure. Throughout most of 
the structure, p is lower than unity, showing again that magnetic 
support is ubiquitous in the system. As examined in the previous 
section, for AD the pile-up of the field occurs outside the core 
and is much less extended (see Fig. 4, at r ~ 30 au), and is con¬ 
trolled by the physical, rather than numerical, resistivity. Further¬ 
more, it does not lead to the growth of the pseudo-disk^. Fig. 17 
compares the density/velocity and plasma p maps, viewed from 


^ Growth of the pseudo-disk due to pile-up of the field has been studied 
in Hennebelle & Fromang (2008) in the appendix with a model based on 
quasi-equilibrium growth of a magnetised self-gravitating rotationally 
dominated structure in iMHD. 



Fig. 15. Evolution of the disk mass (solid line) and the peak magnetic 
field value (dashed line) for the yu = 5, aligned case, in iMHD. The 
four red dots marked 1 to 4 indicate the times at which the snapshots 
in Fig. 16 were taken. We recall that the formation of the first Larson 
core occurs at t = 24300 years (see Table 2). The number 1 indicates 
the beginning of the release of magnetic flux from inside the core. 


above, in the iMHD (top) and AD (bottom) cases. We clearly see 
the contrast between the highly magnetised {p <^ \) disk-like 
structure that arises in iMHD and the rotationally supported disk 
that develops spiral arms with A ^ 1 when ambipolar diffusion 
is taken into account. 

5.3. Numerical convergence 

Numerical convergence is always a problem in numerical sim¬ 
ulations. We carried out a resolution study in the aligned case 
and, while good agreement was found between the various AD 
runs, convergence in iMHD was not satisfactory. Numerical dif¬ 
fusion effects, as well as the choice of Riemann solver to com¬ 
pute fluxes at the cell interfaces, prevent a clear identification 
of the source of this problem. This issue has been highlighted 
in Commer 9 on et al. (2010); Hennebelle & Fromang (2008) for 
the influence of different solvers, and explored in detail by Li 
et al. (2014b), to which we refer to for further information and 
discussions. As part of our resolution study, we performed new 
runs of our benchmark aligned p = 5 case (hereafter denomi¬ 
nated case-0), with a mesh resolution increased by a factor 2 
(16 cells per Jeans length) and a resolution decreased to six cells 
per Jeans length (this latter is only applied for densities above 
10“^^ g cm“^, while keeping the original Jeans length sampling 
elsewhere). 

By comparing the high-resolution iMHD calculation to 
case-0, we found a good agreement until t - 25 kyears (the 
disk mass is slightly lower by a factor 1.5 at this stage, but this 
probably stems from the uncertainties of the criteria used to de¬ 
fine the disk when no clearly identifiable rotationally supported 
structures are present; see Sect. 4), at which point magnetic field 
oscillations at the centre of the core start to appear due to the 
flux freezing condition (earlier in the high-resolution run). This 
is reassuring in the sense that it implies that the resolution we 
have used in case-0 is high enough for numerical resistivity ef¬ 
fects to be negligible, which agrees with the absence of a strong 
field accumulation at the boundary of the core. The higher reso¬ 
lution, however, further diminishes the numerical resistivity, and 
thus flux freezing and flux accumulation are increased, facilitat- 


Article number, page 14 of 20 





























J. Masson et al.: Ambipolar diffusion in low-mass star formation. L 


150 

100 

50 


-50 

-100 

-150 

0 20 40 60 80100120140160180 20 40 60 80100120140160180 20 40 60 801001201401601800 20 40 60 80100120140160180 
R (a.U.) time = 26.340 Kyears time = 26.961 Kyears time = 28.406 Kyears time = 31.582 Kyears 




■o 


(Q 

O 


(a) Evolution of the density/velocity map for the four times labelled 1-4 
in Fig. 15. 
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(b) The three first columns represent the relative importance of each 
component of the magnetic field. The colour scale shows the value of 
each component as follows: black for high values (> 1), white for low 
values (< 10"^ ^); blue and red (in-between > 10“^). The rightmost col¬ 
umn displays the plasma parameter yS = T’/Pmag- 


Fig. 16. Snapshots for different times of interest in the iMHD ji - 5 
aligned case, as labelled Fig. 15. 


ing the development of instabilities, which now appear earlier in 
the simulation. 

In the low-resolution run, the numerical resistivity is in¬ 
creased, and we observe the appearance of a diffusion plateau, 
similar to the one in AD runs, but developing at a slightly higher 
value (P ^ 1 G). In this case, the disk mass is overestimated 
but the symmetry is preserved longer than in the high-resolution 
runs. We observe transient ejections similar to those produced 
by the interchange instability, however, that are due to a strong 
pile-up of the magnetic field close to the core. These ejections 
ultimately lead to a similar quasi-steady final state. 

A visual comparison of the different resolutions (6, 8, and 
16 cells per Jeans length) at two different times is shown in 
Fig. 18. As seen, the resolution strongly affects the simulation 



time = 26.524 Kyears time = 26.524 Kyears 

(a) iMHD 



time = 26.511 Kyears time = 26.511 Kyears 

(b) niMHD 

Fig. 17. yu = 5, aligned case. View of the disk plane in a cylindrical 
volume of height h=20 au. Left and right: density/velocity field and 
plasma p, respectively, with the velocity field. 
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Fig. 18. iMHD, yu = 5, aligned case. View of the disk plane in a cylin¬ 
drical control volume of height h=20 au. First row: t ~ 25.6 k years. 
Second row: t ~ 27 k years. From left to right, the same simulation with 
increased resolution: 6, 8, and 16 cells per Jeans length. 


output when numerical resistivity is not negligible. This is of 
prime importance in turbulent simulations, where it is never clear 
that the turbulent reconnection accurately describes the sub-grid 
unresolved physics. Another parameter that leads to numerical 
diffusivity is the numerical method chosen to solve the induction 
equation. More details on this issue are given in Appendix B. 
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Time (kyears) after first core formation 


Fig. 19. niMHD runs. Comparison of the disk mass in case-0 for dif¬ 
ferent numerical resolutions using the HLLD solver (black) and HLL 
solver (red). 


The numerical convergence is much better for niMHD. 
Fig. 19 shows the disk mass evolution using different mesh reso¬ 
lutions and Riemann solvers. The long-term evolutions are simi¬ 
lar, with both the same global variations and final disk mass. The 
differences also remain small when changing the solver from 
HLLD to HLL, as long as the resolution is > 8 cells per Jeans 
length. We used the HLLD solver for all our simulations to en¬ 
sure that ambipolar diffusion dominates the numerical one. 

The conclusion of our convergence study is as follows. While 
in ideal MHD the isothermal first phase of the collapse is accu¬ 
rately described, with the formation of the first core and the early 
growth of the magnetic bubble/tower, the later evolution strongly 
depends on the numerical resolution. The key issue is the sud¬ 
den release of the magnetic energy accumulated in the core that 
leads to oscillations in the peak magnetic field value (see Fig. 15 
dashed blue line). The release of this energy, depending on the 
numerical tools, ultimately leads to either a quasi-steady state 
in which a massive disk-like rotationally and magnetically sup¬ 
ported structure forms (see Fig. 17 in the top panels) or, for 
HLLD, to the appearance of a high-velocity jet and high Alfven 
speeds that are due to depleted cavities above and below the core 
(see Appendix B for a discussion of the Riemann solver). The 
long-term evolution of these structures, which directly inherit 
their properties from the angular momentum and magnetic flux 
repartition close to the core, is thus of highly questionable valid¬ 
ity. In misaligned configurations, the pseudo-disk is thicker and 
the spurious effects due to very high gradients are less acute (see 
e.g. Joos et al. 2012). The interchange instability is less visible 
than in Fig. 18, but the flux release and the growth of a rota¬ 
tionally supported structure (as in Fig. 15) is still present in the 
simulation. 

In niMHD, in contrast, neither the resolution (as long as it is 
good enough) nor the choice of the solver have significant effects 
on the final outputs. This brings confidence to our global study 
and reinforces our general conclusion that ambipolar diffusion 
in the protostellar first collapse plays a major role in regulating 
angular momentum transport and magnetic flux diffusion, not to 
mention preventing spurious instabilities found in ideal MHD. 


6. Comparison with previous work 

In this section, we compare our results with those of four pre¬ 
vious studies that are representative of the exploration of star 
formation in non-ideal magnetohydrodynamics. 

Krasnopolsky et al. (2010) first explored the impact of en¬ 
hanced resistivities on disk formation. Their study was motivated 
by the fact that in ideal MHD, magnetic fields prevent the forma¬ 
tion of large Keplerian disks, which were ubiquitous in hydrody- 
namical simulations. Resistive MHD was identified as a physical 
process able to limit this effect. Their setup uses an isothermal 
equation of state, does not account for self-gravity, and assumes 
axisymmetry. Their main results, compared to ours, are the fol¬ 
lowing: 

- in their Figs. 5-8, the pinching of the field lines in the equator 
and therefore the split monopole configuration is released as 
the resistivity is increased. 

- they found that a resistivity >10^^ cm^ s“^ is needed to form 
a substantial disk. 

While they did not directly link the release of the split-monopole 
configuration to the formation of rotationally supported disks, 
they concluded that increasing the resistivity helps both to form 
disks and to reduce numerical artefacts through physical rather 
than numerical dissipation, and it also changes the topology of 
the field that is prone to reconnect close to the star where the 
pinching was strong. We agree with these conclusions and note 
that the resistivity above which they formed disks corresponds 
to or is slightly stronger than the typical values of ambipolar 
resistivity we used (see Marchand et al. 2015) that lead to the 
formation of disks. However, these authors assumed a constant 
enhanced resistivity in the induction equation in the Laplace op¬ 
erator. Consequently, they could not grasp any non-linear effect, 
which, as we showed, can lead to a saturation of the magnetic 
field and further facilitate the formation of rotationally supported 
structures. 

Li et al. (2014b) continued the work started by Krasnopol¬ 
sky et al. (2010) and studied the mechanisms of disk formation. 
While they focused on iMHD, they raised many of the questions 
we developed in Sect. 5. In particular, they insisted on the role of 
the warping and the reduced reconnection close to the protostar, 
and they emphasized that a correct treatment of magnetic flux ac¬ 
cretion on the protostar (especially if a sink particle is used) and 
reconnection at the boundary of the disk/pseudo-disk are manda¬ 
tory and are lacking in many previous studies. The problem of 
turbulent reconnection in iMHD was also discussed and placed 
in perspective. We agree with their findings and insist on the fact 
that numerical questions regarding flux conservation and recon¬ 
nection are crucial for a proper study of disk formation. 

Tsukamoto et al. (2015) presented the first second core cal¬ 
culations with both Ohmic dissipation and ambipolar diffusion. 
They performed three-dimensional simulations using an SPH 
code with a Godunov module and divergence cleaning methods 
for the induction equation. Their main conclusions are the fol¬ 
lowing: 

- 1) a value ySpiasma ^ lO'^ in a flattened first core. 

- 2) formation of a circumstellar disk (with a radius of 1 au) at 
the formation of the protostar. 

- 3) occurrence of a saturation plateau at 5 ~ a few 10“^ G in 

the evolution of the magnetic field at densities 10‘ < p < 

lO"!"! g cm“5 followed by a sheet-like collapsing phase (in 
which 5 oc pb^) until the formation of the first core. 

- 4) they did not observe disks around the first Larson core. 
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Their values for ySpiasma (point 1) are similar to ours (see Ta¬ 
ble 2). So far, we cannot comment on point 2 since we did not 
address the second core formation in this paper. When studying 
the evolution of the central magnetic field as a function of den¬ 
sity, they only plotted the data at the centre of the system (while 
we showed every cell in the simulation), and they thus failed to 
report a saturation of the magnetic field that is seen in our Fig. 1 . 
Through private communication, we have obtained the confirma¬ 
tion that they indeed observed a diffusion plateau that was due 
to ambipolar diffusion, but its origin or consequences are not 
discussed in their work. We note that they assumed a sheet-like 
accretion through the disk, whereas in our simulations accretion 
onto the first core mostly occurs in channels driven by the hour¬ 
glass configuration. There is, however, no real evidence for their 
assumption. Finally, the lifetime of the first core in their simula¬ 
tions is comparable to the time of integration after the first core 
formation in our runs (of the order of a thousand years for yu = 4). 
However, they did not form any mid-scale (a dozen to a few 
dozen au) rotationally supported structure around the first core. 
Instead, in their Ohmic and ambipolar-rOhmic calculations, the 
flattened first core itself evolves into a rotationally dominated 
structure that they called a disk. These differences in the dynam¬ 
ics of the collapse highlight the difficulty of conducting reliable 
numerical calculations during prestellar core collapse, which in¬ 
volves a wide range of temporal and spatial scales, and the ne¬ 
cessity of including all relevant physical processes and of using 
reliable numerical schemes. These studies, and the present ones, 
open the route to such explorations. 

Tomida et al. (2015) performed a study similar to ours, fo¬ 
cusing on the first core formation and the surrounding structures, 
with the addition of radiation transfer instead of a barotropic 
equation of state. They used only one value of magnetisation 
that is close to our low-magnetisation case. Their model with 
ambipolar diffusion also included Ohmic dissipation, and they 
did not study the effects of ambipolar diffusion alone. We agree 
with the following of their results: 

- a flattening of the first core, especially in their OA 
(Ohmic -I- ambipolar diffusion) run. 

- depleted polar cavities where the ambipolar action is domi¬ 
nant. 

- in the simulation labelled I (ideal MHD), they found counter¬ 
rotation of the core (as visible in their Fig. 4), as in the 
present Fig. 14, but explained it by the interchange insta¬ 
bility mixing the physical quantities. While we also found 
counter-rotation in iMHD, it occurs before the interchange 
instability develops. Following Galli et al. (2006), we explain 
this counter-rotation by a strong pile-up of the magnetic field 
yielding enough magnetic braking (or torque) to completely 
stop the rotation of the core and spin it backwards. The gen¬ 
eral conclusion, however, remains similar: iMHD models are 
unable to accurately describe angular momentum transport 
after the instability has developed. 

- an increase of the mass-to-fiux ratio up to or above a fac¬ 
tor 10 due to Ohmic and Ohmic-rambipolar diffusion action, 
especially at the first core scales. 

“ no splitting of the field lines in the midplane. 

- the use of a genuine physical dissipation scale helps numer¬ 
ical convergence and thus assesses the reliability of the re¬ 
sults. 

Some of their results are more difficult to discuss. For instance, 
their explanation for a higher first core mass in the resistive cases 
where the magnetic pressure and the thermal and rotational sup¬ 
ports are higher. In our experience, the interplay between mag¬ 
netic braking (yielding stronger or weaker rotational support) 


and removal of magnetic flux (yielding stronger or weaker mag¬ 
netic support) strongly affects the mass of the core. Although it 
is possible that in their set of simulation, the iMHD models lead 
to less massive first cores, in our runs we found that the first core 
mass is higher in the iMHD framework for ju = 2, but somewhat 
similar to or lower than the niMHD value for yu = 5 (see Table 2). 
The physical explanation remains the same: less support in total, 
but it is not straightforward to predict the results over a broad 
range of initial conditions. Although their discussion of the disk 
is not expanded enough to allow detailed comparison with our 
results, the main result is similar: non-ideal MHD enables disk 
formation with no need to artificially relax the flux-freezing con¬ 
straint or to invoke enhanced resistivities. 

7. Conclusions 

We have thoroughly explored the role of ambipolar diffusion in 
magneto-hydrodynamic collapse calculations of a dense molec¬ 
ular cloud core in the context of star formation. Our setup in¬ 
volved a magnetised core of uniform density in solid-body ro¬ 
tation, which collapsed under its own gravity. The gas had a 
barotropic equation of state to mimic radiative transfer, and the 
magnetic resistivities entering the calculations of the non-ideal 
MHD terms were calculated using a reduced chemical network. 
We performed eight simulations, with two different magnetisa¬ 
tions (fi = 5 and 2) and two tilting values of the rotation axis 
(0° and 40°) with respect to the magnetic field direction, both 
in ideal MHD and in non-ideal MHD with ambipolar diffusion. 
We paid particular attention to the problems of magnetic flux 
conservation, magnetic braking, and flux release due to diffusion 
processes. Our main findings are summarised as follows. 

• Ambipolar diffusion creates a magnetic diffusion barrier at 
about the time the first Larson core forms, preventing the 
magnetic field from being amplified above 0.1 G. Flux freez¬ 
ing, however, still holds during the initial stages of the col¬ 
lapse, when resistivities are low. The mass and radius of the 
first Larson core remain rather similar between iMHD and 
niMHD models. 

• The magnitude of B at the diffusion plateau can be esti¬ 
mated by simple order-of-magnitude arguments. This satu¬ 
ration value appears to depend very weakly on the initial 
cloud magnetisation, suggesting a convergence of the final 
state once the initial conditions have been "forgotten" by the 
system. 

• The occurrence of a diffusion plateau has crucial conse¬ 
quences on magnetic braking processes. Not only does it pre¬ 
vent a catastrophic amplification of the field B, which con¬ 
trols the braking efficiency, but it also reorganises the field 
topology, reducing the pinching of field lines close to the pro- 
tostellar object, which also hinders the braking mechanism. 

• Magnetic flux freezing and magnetic braking play a central 
role in the formation and development of the structures dur¬ 
ing the collapse. While in iMHD strongly amplified fields 
launch powerful magnetically supported outflows, these lat¬ 
ter are much weaker or even disappear in niMHD. 

• Misalignment between the initial rotation axis and the mag¬ 
netic field direction does not appear to affect the niMHD re¬ 
sults, showing that the physical flux dissipation due to am¬ 
bipolar diffusion dominates the effects of initial configura¬ 
tion or numerical diffusion. For iMHD models, the additional 
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mixing due to a tilted configuration also produces a diffusion 
plateau, similar to the ambipolar calculations. 

• The disks that form in iMHD, characterised by strong mag¬ 
netic support and an inflated shape (due to toroidal field lines 
loaded by infalling matter) strongly differ from the fiat disks 
that form in niMHD. Formation and long-term evolution of 
disks in ideal MHD, however, are of dubious validity. Fur¬ 
thermore, the excessive magnetic braking generates unphysi¬ 
cal counter-rotation inside the outflow or the magnetic tower. 
Interchange instabilities develop at the interface between the 
core and the disk, producing a redistribution of the flux and 
displacement of the core and disrupting the top-bottom sym¬ 
metry in the system. Mesh resolution is found to strongly af¬ 
fect the simulation results in iMHD. None of these effects is 
observed in the niMHD simulations if the resolution is high 
enough. 

• Disks with Keplerian velocity profiles around the protostar 
form in all our niMHD simulations for all different magneti¬ 
sations and inclination angles. Their size and mass, however, 
is significantly reduced (by a factor ~ 10) in the more mag¬ 
netised case (fi = 2) because of the increased braking. Such 
a magnetisation value seems to be typical of most molecular 
clouds. Aligned and misaligned initial configurations have 
no consequence on the disk properties in niMHD and yield 
disks with very similar properties. 

Magnetic flux diffusivity due to ambipolar diffusion thus ap¬ 
pears to play a dominant role during the first protostellar collapse 
and the formation of the first Larson core and its surrounding 
disk, and it appears to dominate processes such as magnetic field 
and rotation axis orientations. Flux diffusion during the collapse 
allows the formation of quasi-Keplerian disks, solving the "disk 
formation crisis" found in ideal MHD. The mass, size, and mag¬ 
netic properties (plasma JS) of these disks, however, strongly de¬ 
pend on the initial magnetisation (p parameter) of the cloud, and 
in all cases the disks are significantly smaller and less massive 
than those found in pure hydrodynamics calculations. Charac¬ 
terisation of this plasma p is of major importance to characterise 
viscous transport in disks, notably in the MRI. To complete our 
study, we will explore the effect of turbulence on the first col¬ 
lapse of a magnetised cloud core in a forthcoming paper. 
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Appendix A: Calculation of the saturation value 


In this section, we provide an analytical estimate of the value 
of B at which the magnetic diffusion starts to overcome the dy¬ 
namical MHD effects. We rewrite the initial values for and 
Po in terms of the ratio of thermal to gravitational energy a and 
mass-to-fiux ratio p, according to 



(A.l) 

(A.2) 



Fig. B.2. First core mass evolution for various solvers. Standard resolu¬ 
tion of 8 cells per Jeans length. 


The relationship between the magnetic field and density is as¬ 
sumed to follow p oc We also assume that in the early phases 
of collapse, the typical length scale is L ^ <^sfiree-faii and the ve¬ 
locity is y ^ Cg. We then seek the value of the magnetic field ^sat 
for which 


VxB=.,„ix 




(A.3) 


with Pad - y^^p.p • Using the fact that pi = C ^ (from Shu 
1992), we can then write 

B^ = Cs^free-fallTADP^^^U, (A.4) 

which yields the saturation value for the magnetic field 


^sat ~ UyAD 


3 

InG 


A 




(A.5) 


Appendix B: Influence of the solver 

The choice of the Riemann solver to compute the conservative 
fluxes at the interfaces between domain cells can greatly affect 
the results of a simulation, especially when studying diffusion 
processes. We compare here four test simulations that use dif¬ 
ferent strategies to compute the fluxes. The resolution for each 
run is 8 cells per Jeans length. In Fig. B.l from left to right, the 
solver is increasingly less diffusive. 

In the HLLD-l-switch simulation, we used the HLLD solver 
over most of the computational domain, but switched to a Lax- 
Friedrich solver whenever large discontinuities^ are present in 
the flow variables between neighbouring cells. This allowed us 
to evolve the simulations for long periods of time, up to ~ 1.5 
free-fall time. However, it raises crucial problems on both the 
numerical reconnection in the vicinity of the core and the conse¬ 
quent interchange instability, as well as a violent core displace¬ 
ment that accompanies an unphysical release of magnetic energy 
at the centre. The formation of large inflated disks is also ques¬ 
tionable. 

In a second calculation, we removed the switch to a more dif¬ 
fusive solver, keeping HLLD in the entire domain. In this case, 

^ Discontinuities either in density, or when the wave speed for the in¬ 
termediate states for HLLD are ill defined due to a very small denomi¬ 
nator. 
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Fig. B.l. Density maps of the region around the first core using four different Riemann solvers: (a) Lax-Friedrich, (b) HLL, (c) HLLD-i-Switch, 
(d) HLLD. The velocity vectors are overlaid on density maps with green and blue arrows. 


the discontinuities (mostly magnetic and density related) con¬ 
tinue to sharpen, until the time integration in the simulation is 
frozen due to very high Alfven speeds in the outflows or jets (it 
can reach speeds of the order of > 10^ km s“^). Whether this is a 
correct behaviour is open to debate (see also Li et al. 2014b). We 
reach convergence of results quicker than in the HLLD-h switch 
case, since we basically prevent numerical diffusion, provided 
we have enough resolution at a given scale and density. On 
the other hand, there are physical diffusive processes at play in 
star formation, either from microphysics (e.g. ambipolar diffu¬ 
sion) or from turbulent reconnection. Therefore, a description 
that avoids any diffusion lacks physical mechanisms. Finally, we 
performed simulations using only HLL that takes into account 
only three waves instead of five with HLLD. This case fits in 
between the two previous ones; convergence is reached with 16 
cells per Jeans length, but we do witness diffusion in the out¬ 
flow and the vicinity or interior of the Larson core. In this case, 
Alfven speeds remain below 10^ km s“^ For simulations of the 
second collapse and formation of the second core, for which time 
integration is critical, our former tests suggest that using HLL 
with high enough resolution, at least in ideal MHD, is better than 
using HLLD (with or without the switch) because it enables both 
numerical convergence and the probing of longer timescales. For 
the sake of comparison, we performed one test case using the 
Lax-Friedrich solver alone. The results are similar to the HLL 
case. 


Additionally, Fig. B.2 shows the first core mass evolution for 
each simulation, highlighting significant differences in an out¬ 
come as well defined as the flrst core mass. Until the time when 
the least diffusive solver (HLLD) does not allow continuing the 
simulation anymore because the time step is too small, the results 
are similar. However, after this point, depending on the choice of 
solver, results can significantly differ both in the qualitative pic¬ 
ture (as seen in Fig. B.l) or for a precise outcome (Fig. B.2). 
Therefore, depending on the timescale of interest and the pre¬ 
cise point of study, a careful choice of the numerical method has 
to be made. 


Appendix C: Disk veiocity profiie using ideai MHD 

In this appendix, we caution about the definition of a Keplerian 
disk and the estimate of the central core mass obtained when 
fitting a priori the disk radial velocity distribution with a Kep¬ 
lerian profile (cx in particular in structures whose growth 

is magnetically driven, as is the case for pseudo-disks in ideal 
MHD. The disk criteria we used, based on Joos et al. (2012), 
define a rotationally dominated structure. We then proceeded by 
studying the radial dependence of the toroidal velocity, and com¬ 
pared it to the classical Keplerian velocity profiles. The velocity 
profile, the plasma (3, and the aspect ratio allowed us to estimate 
how similar the disks we found are to a Keplerian disk. As em¬ 
phasised below, we And that these disks can depart significantly 
from the classical Keplerian picture. The relatively high mass of 
the surroundings of the flrst core and the magnetic fields cause 
the usual Keplerian model to fail. This issue is of prime impor¬ 
tance when trying to characterise the properties of the disk struc¬ 
tures observed around Class-0 objects (see e.g. Tobin et al. 2013, 
2015). 

An example is shown in Fig. C.l. Whereas a Keplerian pro¬ 
file (grey dashed line) can be roughly fitted to the toroidal ve¬ 
locity held (blue dots), the best-fitting estimate (dashed red line) 
gives a significantly different value, cx yielding a central 

mass different from the one inferred from the model of a Ke¬ 
plerian velocity profile around a central point mass. This result 
stresses the fact that fitting a scattered velocity profile, which can 
depart significantly from a Keplerian one, by such a value can be 
very uncertain, casting doubts on estimates of central masses ob¬ 
tained with Keplerian formulae. 
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